{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecasting in statsmodels\n",
    "\n",
    "This notebook describes forecasting using time series models in statsmodels.\n",
    "\n",
    "**Note**: this notebook applies only to the state space model classes, which are:\n",
    "\n",
    "- `sm.tsa.SARIMAX`\n",
    "- `sm.tsa.UnobservedComponents`\n",
    "- `sm.tsa.VARMAX`\n",
    "- `sm.tsa.DynamicFactor`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "macrodata = sm.datasets.macrodata.load_pandas().data\n",
    "macrodata.index = pd.period_range('1959Q1', '2009Q3', freq='Q')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic example\n",
    "\n",
    "A simple example is to use an AR(1) model to forecast inflation. Before forecasting, let's take a look at the series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endog = macrodata['infl']\n",
    "endog.plot(figsize=(15, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constructing and estimating the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to formulate the econometric model that we want to use for forecasting. In this case, we will use an AR(1) model via the `SARIMAX` class in statsmodels.\n",
    "\n",
    "After constructing the model, we need to estimate its parameters. This is done using the `fit` method. The `summary` method produces several convenient tables showing the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the model\n",
    "mod = sm.tsa.SARIMAX(endog, order=(1, 0, 0), trend='c')\n",
    "# Estimate the parameters\n",
    "res = mod.fit()\n",
    "\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Out-of-sample forecasts are produced using the `forecast` or `get_forecast` methods from the results object.\n",
    "\n",
    "The `forecast` method gives only point forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The default is to get a one-step-ahead forecast:\n",
    "print(res.forecast())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `get_forecast` method is more general, and also allows constructing confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we construct a more complete results object.\n",
    "fcast_res1 = res.get_forecast()\n",
    "\n",
    "# Most results are collected in the `summary_frame` attribute.\n",
    "# Here we specify that we want a confidence level of 90%\n",
    "print(fcast_res1.summary_frame(alpha=0.10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default confidence level is 95%, but this can be controlled by setting the `alpha` parameter, where the confidence level is defined as $(1 - \\alpha) \\times 100\\%$. In the example above, we specified a confidence level of 90%, using `alpha=0.10`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specifying the number of forecasts\n",
    "\n",
    "Both of the functions `forecast` and `get_forecast` accept a single argument indicating how many forecasting steps are desired. One option for this argument is always to provide an integer describing the number of steps ahead you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.forecast(steps=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcast_res2 = res.get_forecast(steps=2)\n",
    "# Note: since we did not specify the alpha parameter, the\n",
    "# confidence level is at the default, 95%\n",
    "print(fcast_res2.summary_frame())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, **if your data included a Pandas index with a defined frequency** (see the section at the end on Indexes for more information), then you can alternatively specify the date through which you want forecasts to be produced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.forecast('2010Q2'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcast_res3 = res.get_forecast('2010Q2')\n",
    "print(fcast_res3.summary_frame())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the data, forecasts, and confidence intervals\n",
    "\n",
    "Often it is useful to plot the data, the forecasts, and the confidence intervals. There are many ways to do this, but here's one example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(15, 5))\n",
    "\n",
    "# Plot the data (here we are subsetting it to get a better look at the forecasts)\n",
    "endog.loc['1999':].plot(ax=ax)\n",
    "\n",
    "# Construct the forecasts\n",
    "fcast = res.get_forecast('2011Q4').summary_frame()\n",
    "fcast['mean'].plot(ax=ax, style='k--')\n",
    "ax.fill_between(fcast.index, fcast['mean_ci_lower'], fcast['mean_ci_upper'], color='k', alpha=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note on what to expect from forecasts\n",
    "\n",
    "The forecast above may not look very impressive, as it is almost a straight line. This is because this is a very simple, univariate forecasting model. Nonetheless, keep in mind that these simple forecasting models can be extremely competitive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction vs Forecasting\n",
    "\n",
    "The results objects also contain two methods that all for both in-sample fitted values and out-of-sample forecasting. They are `predict` and `get_prediction`. The `predict` method only returns point predictions (similar to `forecast`), while the `get_prediction` method also returns additional results (similar to `get_forecast`).\n",
    "\n",
    "In general, if your interest is out-of-sample forecasting, it is easier to stick to the `forecast` and `get_forecast` methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross validation\n",
    "\n",
    "**Note**: some of the functions used in this section were first introduced in statsmodels v0.11.0.\n",
    "\n",
    "A common use case is to cross-validate forecasting methods by performing h-step-ahead forecasts recursively using the following process:\n",
    "\n",
    "1. Fit model parameters on a training sample\n",
    "2. Produce h-step-ahead forecasts from the end of that sample\n",
    "3. Compare forecasts against test dataset to compute error rate\n",
    "4. Expand the sample to include the next observation, and repeat\n",
    "\n",
    "Economists sometimes call this a pseudo-out-of-sample forecast evaluation exercise, or time-series cross-validation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will conduct a very simple exercise of this sort using the inflation dataset above. The full dataset contains 203 observations, and for expositional purposes we'll use the first 80% as our training sample and only consider one-step-ahead forecasts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single iteration of the above procedure looks like the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: fit model parameters w/ training sample\n",
    "training_obs = int(len(endog) * 0.8)\n",
    "\n",
    "training_endog = endog[:training_obs]\n",
    "training_mod = sm.tsa.SARIMAX(\n",
    "    training_endog, order=(1, 0, 0), trend='c')\n",
    "training_res = training_mod.fit()\n",
    "\n",
    "# Print the estimated parameters\n",
    "print(training_res.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: produce one-step-ahead forecasts\n",
    "fcast = training_res.forecast()\n",
    "\n",
    "# Step 3: compute root mean square forecasting error\n",
    "true = endog.reindex(fcast.index)\n",
    "error = true - fcast\n",
    "\n",
    "# Print out the results\n",
    "print(pd.concat([true.rename('true'),\n",
    "                 fcast.rename('forecast'),\n",
    "                 error.rename('error')], axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To add on another observation, we can use the `append` or `extend` results methods. Either method can produce the same forecasts, but they differ in the other results that are available:\n",
    "\n",
    "- `append` is the more complete method. It always stores results for all training observations, and it optionally allows refitting the model parameters given the new observations (note that the default is *not* to refit the parameters).\n",
    "- `extend` is a faster method that may be useful if the training sample is very large. It *only* stores results for the new observations, and it does not allow refitting the model parameters (i.e. you have to use the parameters estimated on the previous sample).\n",
    "\n",
    "If your training sample is relatively small (less than a few thousand observations, for example) or if you want to compute the best possible forecasts, then you should use the `append` method. However, if that method is infeasible (for example, because you have a very large training sample) or if you are okay with slightly suboptimal forecasts (because the parameter estimates will be slightly stale), then you can consider the `extend` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A second iteration, using the `append` method and refitting the parameters, would go as follows (note again that the default for `append` does not refit the parameters, but we have overridden that with the `refit=True` argument):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: append a new observation to the sample and refit the parameters\n",
    "append_res = training_res.append(endog[training_obs:training_obs + 1], refit=True)\n",
    "\n",
    "# Print the re-estimated parameters\n",
    "print(append_res.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that these estimated parameters are slightly different than those we originally estimated. With the new results object, `append_res`, we can compute forecasts starting from one observation further than the previous call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: produce one-step-ahead forecasts\n",
    "fcast = append_res.forecast()\n",
    "\n",
    "# Step 3: compute root mean square forecasting error\n",
    "true = endog.reindex(fcast.index)\n",
    "error = true - fcast\n",
    "\n",
    "# Print out the results\n",
    "print(pd.concat([true.rename('true'),\n",
    "                 fcast.rename('forecast'),\n",
    "                 error.rename('error')], axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting it altogether, we can perform the recursive forecast evaluation exercise as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup forecasts\n",
    "nforecasts = 3\n",
    "forecasts = {}\n",
    "\n",
    "# Get the number of initial training observations\n",
    "nobs = len(endog)\n",
    "n_init_training = int(nobs * 0.8)\n",
    "\n",
    "# Create model for initial training sample, fit parameters\n",
    "init_training_endog = endog.iloc[:n_init_training]\n",
    "mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend='c')\n",
    "res = mod.fit()\n",
    "\n",
    "# Save initial forecast\n",
    "forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Step through the rest of the sample\n",
    "for t in range(n_init_training, nobs):\n",
    "    # Update the results by appending the next observation\n",
    "    updated_endog = endog.iloc[t:t+1]\n",
    "    res = res.append(updated_endog, refit=False)\n",
    "    \n",
    "    # Save the new set of forecasts\n",
    "    forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Combine all forecasts into a dataframe\n",
    "forecasts = pd.concat(forecasts, axis=1)\n",
    "\n",
    "print(forecasts.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a set of three forecasts made at each point in time from 1999Q2 through 2009Q3. We can construct the forecast errors by subtracting each forecast from the actual value of `endog` at that point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the forecast errors\n",
    "forecast_errors = forecasts.apply(lambda column: endog - column).reindex(forecasts.index)\n",
    "\n",
    "print(forecast_errors.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To evaluate our forecasts, we often want to look at a summary value like the root mean square error. Here we can compute that for each horizon by first flattening the forecast errors so that they are indexed by horizon and then computing the root mean square error fore each horizon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reindex the forecasts by horizon rather than by date\n",
    "def flatten(column):\n",
    "    return column.dropna().reset_index(drop=True)\n",
    "\n",
    "flattened = forecast_errors.apply(flatten)\n",
    "flattened.index = (flattened.index + 1).rename('horizon')\n",
    "\n",
    "print(flattened.iloc[:3, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the root mean square error\n",
    "rmse = (flattened**2).mean(axis=1)**0.5\n",
    "\n",
    "print(rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using `extend`\n",
    "\n",
    "We can check that we get similar forecasts if we instead use the `extend` method, but that they are not exactly the same as when we use `append` with the `refit=True` argument. This is because `extend` does not re-estimate the parameters given the new observation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup forecasts\n",
    "nforecasts = 3\n",
    "forecasts = {}\n",
    "\n",
    "# Get the number of initial training observations\n",
    "nobs = len(endog)\n",
    "n_init_training = int(nobs * 0.8)\n",
    "\n",
    "# Create model for initial training sample, fit parameters\n",
    "init_training_endog = endog.iloc[:n_init_training]\n",
    "mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend='c')\n",
    "res = mod.fit()\n",
    "\n",
    "# Save initial forecast\n",
    "forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Step through the rest of the sample\n",
    "for t in range(n_init_training, nobs):\n",
    "    # Update the results by appending the next observation\n",
    "    updated_endog = endog.iloc[t:t+1]\n",
    "    res = res.extend(updated_endog)\n",
    "    \n",
    "    # Save the new set of forecasts\n",
    "    forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)\n",
    "\n",
    "# Combine all forecasts into a dataframe\n",
    "forecasts = pd.concat(forecasts, axis=1)\n",
    "\n",
    "print(forecasts.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the forecast errors\n",
    "forecast_errors = forecasts.apply(lambda column: endog - column).reindex(forecasts.index)\n",
    "\n",
    "print(forecast_errors.iloc[:5, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reindex the forecasts by horizon rather than by date\n",
    "def flatten(column):\n",
    "    return column.dropna().reset_index(drop=True)\n",
    "\n",
    "flattened = forecast_errors.apply(flatten)\n",
    "flattened.index = (flattened.index + 1).rename('horizon')\n",
    "\n",
    "print(flattened.iloc[:3, :5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the root mean square error\n",
    "rmse = (flattened**2).mean(axis=1)**0.5\n",
    "\n",
    "print(rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By not re-estimating the parameters, our forecasts are slightly worse (the root mean square error is higher at each horizon). However, the process is faster, even with only 200 datapoints. Using the `%%timeit` cell magic on the cells above, we found a runtime of 570ms using `extend` versus 1.7s using `append` with `refit=True`. (Note that using `extend` is also faster than using `append` with `refit=False`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexes\n",
    "\n",
    "Throughout this notebook, we have been making use of Pandas date indexes with an associated frequency. As you can see, this index marks our data as at a quarterly frequency, between 1959Q1 and 2009Q3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(endog.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In most cases, if your data has an associated data/time index with a defined frequency (like quarterly, monthly, etc.), then it is best to make sure your data is a Pandas series with the appropriate index. Here are three examples of this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Annual frequency, using a PeriodIndex\n",
    "index = pd.period_range(start='2000', periods=4, freq='A')\n",
    "endog1 = pd.Series([1, 2, 3, 4], index=index)\n",
    "print(endog1.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quarterly frequency, using a DatetimeIndex\n",
    "index = pd.date_range(start='2000', periods=4, freq='QS')\n",
    "endog2 = pd.Series([1, 2, 3, 4], index=index)\n",
    "print(endog2.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monthly frequency, using a DatetimeIndex\n",
    "index = pd.date_range(start='2000', periods=4, freq='M')\n",
    "endog3 = pd.Series([1, 2, 3, 4], index=index)\n",
    "print(endog3.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, if your data has an associated date/time index, it is best to use that even if does not have a defined frequency. An example of that kind of index is as follows - notice that it has `freq=None`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = pd.DatetimeIndex([\n",
    "    '2000-01-01 10:08am', '2000-01-01 11:32am',\n",
    "    '2000-01-01 5:32pm', '2000-01-02 6:15am'])\n",
    "endog4 = pd.Series([0.2, 0.5, -0.1, 0.1], index=index)\n",
    "print(endog4.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can still pass this data to statsmodels' model classes, but you will get the following warning, that no frequency data was found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = sm.tsa.SARIMAX(endog4)\n",
    "res = mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this means is that you cannot specify forecasting steps by dates, and the output of the `forecast` and `get_forecast` methods will not have associated dates. The reason is that without a given frequency, there is no way to determine what date each forecast should be assigned to. In the example above, there is no pattern to the date/time stamps of the index, so there is no way to determine what the next date/time should be (should it be in the morning of 2000-01-02? the afternoon? or maybe not until 2000-01-03?).\n",
    "\n",
    "For example, if we forecast one-step-ahead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.forecast(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index associated with the new forecast is `4`, because if the given data had an integer index, that would be the next value. A warning is given letting the user know that the index is not a date/time index.\n",
    "\n",
    "If we try to specify the steps of the forecast using a date, we will get the following exception:\n",
    "\n",
    "    KeyError: 'The `end` argument could not be matched to a location related to the index of the data.'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Here we'll catch the exception to prevent printing too much of\n",
    "# the exception trace output in this notebook\n",
    "try:\n",
    "    res.forecast('2000-01-03')\n",
    "except KeyError as e:\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ultimately there is nothing wrong with using data that does not have an associated date/time frequency, or even using data that has no index at all, like a Numpy array. However, if you can use a Pandas series with an associated frequency, you'll have more options for specifying your forecasts and get back results with a more useful index."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
